##ROC curves of Table 1, Model 1; Table 1, Model 2; and a model with No Challenger Reputation & No Interaction Term/Treatment
##No single curve appears to strictly dominate the others.
##10 fold cross-validation included below (and not included in our paper)

library(foreign)
library(Zelig)
library(sandwich)
library(MASS)
library(lme4)
library(VGAM)
library(mlogit)

data <- read.dta(file="clare.dta")
head(data)
dim(data)


model1Vars <- c("cwinit", "rivalry_count1", "weakmr_decline9", "weakmr2_decline9", "relcap", "jdem", "pol_rel", "pceyrs", "pceyrs2",  "pceyrs3")

dataNew <- data[, model1Vars]
#head(dataNew)
#nrow(dataNew)
dataNewNoNA <- na.omit(dataNew)
#nrow(dataNewNoNA)

Y<-as.matrix(dataNewNoNA[,1])
#Y<-dataNewNoNA[,1]
X<-as.matrix(cbind(1,dataNewNoNA[,-1]))


table1Model1 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = dataNewNoNA)

table1Model1BETAs<-summary(table1Model1)$coefficients[,1]

fitted.Y <- pnorm(X%*%table1Model1BETAs)

cutoff <- c(seq(0, .1, by=.0001), seq(.1, 1, by=.005))


correct0s <- c()
correct1s <- c()
for (i in 1:length(cutoff)) {
  predictions <- as.numeric(fitted.Y > cutoff[i])
  correct0s[i] <- mean(predictions[Y==0]==0)
  correct1s[i] <- mean(predictions[Y==1]==1)
}

table1Model1correct1s <- correct1s
table1Model1correct0s <- correct0s


##
table1Model2Vars <- c("cwinit", "rivalry_count1", "weakmr_decline9", "nriv_weakdec9mr", "weakmr2_decline9", "relcap", "jdem", "pol_rel", "pceyrs", "pceyrs2",  "pceyrs3")

dataNew <- data[, table1Model2Vars]
dataNewNoNA <- na.omit(dataNew)

Y<-as.matrix(dataNewNoNA[,1])
#Y<-dataNewNoNA[,1]
X<-as.matrix(cbind(1,dataNewNoNA[,-1]))

table1Model2 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = dataNewNoNA)
#coef(table1Model2)


table1Model2BETAs<-summary(table1Model2)$coefficients[,1]

fitted.Y <- pnorm(X%*%table1Model2BETAs)

cutoff <- c(seq(0, .1, by=.0001), seq(.1, 1, by=.005))


correct0s <- c()
correct1s <- c()
for (i in 1:length(cutoff)) {
  predictions <- as.numeric(fitted.Y > cutoff[i])
  correct0s[i] <- mean(predictions[Y==0]==0)
  correct1s[i] <- mean(predictions[Y==1]==1)
}

table1Model2correct1s <- correct1s
table1Model2correct0s <- correct0s



##
noTreatmentModelVars <- c("cwinit", "weakmr2_decline9", "relcap", "jdem", "pol_rel", "pceyrs", "pceyrs2",  "pceyrs3")

dataNew <- data[, noTreatmentModelVars]
dataNewNoNA <- na.omit(dataNew)

Y<-as.matrix(dataNewNoNA[,1])
#Y<-dataNewNoNA[,1]
X<-as.matrix(cbind(1,dataNewNoNA[,-1]))

noTreatmentModel <- zelig(cwinit ~ weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = dataNewNoNA)
#coef(table1Model2)


noTreatmentModelBETAs<-summary(noTreatmentModel)$coefficients[,1]

fitted.Y <- pnorm(X%*%noTreatmentModelBETAs)

cutoff <- c(seq(0, .1, by=.0001), seq(.1, 1, by=.005))


correct0s <- c()
correct1s <- c()
for (i in 1:length(cutoff)) {
  predictions <- as.numeric(fitted.Y > cutoff[i])
  correct0s[i] <- mean(predictions[Y==0]==0)
  correct1s[i] <- mean(predictions[Y==1]==1)
}

noTreatmentModelcorrect1s <- correct1s
noTreatmentModelcorrect0s <- correct0s





#pdf(file="dropTreatmentROC.pdf") 
plot(table1Model1correct1s, table1Model1correct0s, type="l", xlim=c(0,1), ylim=c(0,1), 
     main="ROC Curves",
	 xlab="Proportion of 1's correctly classified",
	 ylab="Proportion of 0's correctly classified")
abline(a=1, b=-1, lty=2)

lines(x = table1Model2correct1s, y = table1Model2correct0s, col="red")
lines(x = noTreatmentModelcorrect1s, y = noTreatmentModelcorrect0s, col="blue")
legend("topright", c("Table 1, Model 1 (No Interaction Term)","Table 1, Model 2 (Interaction Term)", "No Treatment-Related Terms"), bty = "n",
col=c("black", "red", "blue"), lwd=c(1, 1, 1), lty=c(1, 1, 1))
#dev.off()








######################################################VALIDATION CODE (follows code provided in Gov2001 Section):
data <- read.dta(file="clare.dta")
head(data)
dim(data)


model1Vars <- c("cwinit", "rivalry_count1", "weakmr_decline9", "weakmr2_decline9", "relcap", "jdem", "pol_rel", "pceyrs", "pceyrs2",  "pceyrs3")

dataNew <- data[, model1Vars]
#head(dataNew)
#nrow(dataNew)
dataNewNoNA <- na.omit(dataNew)
#nrow(dataNewNoNA)

Y<-as.matrix(dataNewNoNA[,1])
#Y<-dataNewNoNA[,1]
X<-as.matrix(cbind(1,dataNewNoNA[,-1]))


#10-fold validation:

set.seed(1492)
N <- length(Y) 
size <- floor(N/10) 
results <- list()

for (i in 1:10) {
  if(i==1) training <- (size+1):length(Y)
  if(i>1 && i <10) training <- (1:length(Y))[-c(((i-1)*size+1):((i)*size))] 
  if(i==10) training <- 1:(9*size)
  
  X.train <- X[training,]
  Y.train <- Y[training]
  X.test <- X[-training,]
  #names(Y[1]) <- "cwinit" 
  dataNewNoNA2 <- data.frame(cbind(Y.train, X.train))
  #head(dataNewNoNA2)
  #names(dataNewNoNA2[1]) <- "cwinit"
  modelOutput <- zelig(Y.train ~ rivalry_count1 + weakmr_decline9 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = dataNewNoNA2)

  BETAs<-summary(modelOutput)$coefficients[,1]

  fitted.Y <- pnorm(X.test%*%BETAs) 
  		  
  results[[i]] <- fitted.Y	 
}

fitted.Y <- unlist(results)
#par(mfrow=c(1,2))


cutoff <- c(seq(0, .1, by=.0001), seq(.1, 1, by=.005))


correct0s <- c()
correct1s <- c()
for (i in 1:length(cutoff)) {
  predictions <- as.numeric(fitted.Y > cutoff[i])
  correct0s[i] <- mean(predictions[Y==0]==0)
  correct1s[i] <- mean(predictions[Y==1]==1)
}

table1Model1correct1s <- correct1s
table1Model1correct0s <- correct0s


##OUT OF SAMPLE
table1Model2Vars <- c("cwinit", "rivalry_count1", "weakmr_decline9", "nriv_weakdec9mr", "weakmr2_decline9", "relcap", "jdem", "pol_rel", "pceyrs", "pceyrs2",  "pceyrs3")

dataNew <- data[, table1Model2Vars]
dataNewNoNA <- na.omit(dataNew)

Y<-as.matrix(dataNewNoNA[,1])
#Y<-dataNewNoNA[,1]
X<-as.matrix(cbind(1,dataNewNoNA[,-1]))


#10-fold validation:

set.seed(1492)
N <- length(Y) 
size <- floor(N/10) 
results <- list()

for (i in 1:10) {
  if(i==1) training <- (size+1):length(Y)
  if(i>1 && i <10) training <- (1:length(Y))[-c(((i-1)*size+1):((i)*size))] 
  if(i==10) training <- 1:(9*size)
  
  X.train <- X[training,]
  Y.train <- Y[training]
  X.test <- X[-training,]
  #names(Y[1]) <- "cwinit" 
  dataNewNoNA2 <- data.frame(cbind(Y.train, X.train))
  #head(dataNewNoNA2)
  #names(dataNewNoNA2[1]) <- "cwinit"
  modelOutput <- zelig(Y.train ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = dataNewNoNA2)

  BETAs<-summary(modelOutput)$coefficients[,1]

  fitted.Y <- pnorm(X.test%*%BETAs) 
  	  
  results[[i]] <- fitted.Y	 
}

fitted.Y <- unlist(results)


cutoff <- c(seq(0, .1, by=.0001), seq(.1, 1, by=.005))


correct0s <- c()
correct1s <- c()
for (i in 1:length(cutoff)) {
  predictions <- as.numeric(fitted.Y > cutoff[i])
  correct0s[i] <- mean(predictions[Y==0]==0)
  correct1s[i] <- mean(predictions[Y==1]==1)
}

table1Model2correct1s <- correct1s
table1Model2correct0s <- correct0s



##OUT OF SAMPLE
noTreatmentModelVars <- c("cwinit", "weakmr2_decline9", "relcap", "jdem", "pol_rel", "pceyrs", "pceyrs2",  "pceyrs3")

dataNew <- data[, noTreatmentModelVars]
dataNewNoNA <- na.omit(dataNew)

Y<-as.matrix(dataNewNoNA[,1])
#Y<-dataNewNoNA[,1]
X<-as.matrix(cbind(1,dataNewNoNA[,-1]))

#10-fold validation:

set.seed(1492)
N <- length(Y) 
size <- floor(N/10) 
results <- list()

for (i in 1:10) {
  if(i==1) training <- (size+1):length(Y)
  if(i>1 && i <10) training <- (1:length(Y))[-c(((i-1)*size+1):((i)*size))] 
  if(i==10) training <- 1:(9*size)
  
  X.train <- X[training,]
  Y.train <- Y[training]
  X.test <- X[-training,]
  #names(Y[1]) <- "cwinit" 
  dataNewNoNA2 <- data.frame(cbind(Y.train, X.train))
  #head(dataNewNoNA2)
  #names(dataNewNoNA2[1]) <- "cwinit"
  modelOutput <- zelig(Y.train ~ weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = dataNewNoNA2)

  BETAs<-summary(modelOutput)$coefficients[,1]

  fitted.Y <- pnorm(X.test%*%BETAs) 
  			  
  results[[i]] <- fitted.Y	 
}

fitted.Y <- unlist(results)


cutoff <- c(seq(0, .1, by=.0001), seq(.1, 1, by=.005))


correct0s <- c()
correct1s <- c()
for (i in 1:length(cutoff)) {
  predictions <- as.numeric(fitted.Y > cutoff[i])
  correct0s[i] <- mean(predictions[Y==0]==0)
  correct1s[i] <- mean(predictions[Y==1]==1)
}

noTreatmentModelcorrect1s <- correct1s
noTreatmentModelcorrect0s <- correct0s





#pdf(file="dropTreatmentROC.pdf") 
plot(table1Model1correct1s, table1Model1correct0s, type="l", xlim=c(0,1), ylim=c(0,1), 
     main="ROC Curves (10-Fold Validation)",
	 xlab="Proportion of 1's correctly classified",
	 ylab="Proportion of 0's correctly classified")
abline(a=1, b=-1, lty=2)

lines(x = table1Model2correct1s, y = table1Model2correct0s, col="red")
lines(x = noTreatmentModelcorrect1s, y = noTreatmentModelcorrect0s, col="blue")
legend("topright", c("Table 1, Model 1 (No Interaction Term)","Table 1, Model 2 (Interaction Term)", "No Treatment-Related Terms"), bty = "n",
col=c("black", "red", "blue"), lwd=c(1, 1, 1), lty=c(1, 1, 1))
#dev.off()




